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Quantum optimal control theory (QOCT) aims at finding an external field that drives a quantum 
system in such a way that optimally achieves some predefined target. In practice this normally 
means optimizing the value of some observable, a so called merit function. In consequence, a key 
part of the theory is a set of equations, which provides the gradient of the merit function with 
respect to parameters that control the shape of the driving field. We show that these equations 
can be straightforwardly derived using the standard linear response theory, only requiring a minor 
generalization - the unperturbed Hamiltonian is allowed to be time-dependent. As a result, the 
aforementioned gradients are identified with certain response functions. This identification leads to a 
natural reformulation of QOCT in term of the Keldysh contour formalism of the quantum many-body 
theory. In particular, the gradients of the merit function can be calculated using the diagrammatic 
technique for non-equilibrium Green's functions, which should be helpful in the application of QOCT 
to computationally difficult many-electron problems. 

PACS numbers: 02.30.Yy, 71.10.-w, 32.80.Qk 



I. INTRODUCTION 

Quantum Optimal Control Theory (QOCT) 0, i| is 
concerned with finding a time-dependent external field 
that drives a given quantum system to optimally achieve 
some predefined target, that depends on the manner in 
which the system evolves Q. For example, a target can 
be the population of some excited state at the final time 
of the propagation - but many other options are pos- 
sible. The theory can be regarded as a branch of the 
classical control theories developed mostly in the fields 
of mathematics and engineering 0, The quantum 
discipline was born in the late 80s as the most 

complete theoretical framework capable of addressing the 
nascent experimental field of quantum control (or coher- 
ent control) 0. The range of applications of QOCT is 
growing very fast, thanks to the progress in the ultrafast 
laser pulse generation and pulse shaping techniques [l3 , 
as well as to the development of adaptive feedback con- 
trol schemes [ll|, ES] ■ Typical examples of applications 
are the control of the population of excited states in 
molecules 12|, optimization of high-harmonic genera- 
tion [isj . optimization of selective photo-dissociation of 
molecules Il4i , optimization of multi-photon ionization 
of atoms [l5|, enhancement of electron transfer in dye- 
sensitized solar cells [1^1, etc. 

At the formal level, the central problem of QOCT is 
to maximize an expectation value of some operator, usu- 
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ally known as a metrit (or target) function, whose input 
is the external field that needs to be optimally shaped. 
The field is normally parameterized either by a decrete 
set of real-valued "control" parameters, or, in a more 
general setting, by continuous functions of time. In the 
latter case, one usually speaks of target functionals. In 
most cases, the optimization algorithm will require both 
the computation of the merit function and of its gra- 
dient with respect to control parameters. Therefore an 
expression and computational strategy for this gradient 
constitutes one of the most important parts of QOCT. 

The usual derivation of expressions for the gradient 
of the merit function proceeds via the definition of a 
Lagrangian functional, and of a "Lagrange multiplier" 
wave function (see, for example, Refs. 0, [l3)- It leads 
to an expression for the gradient that involves the for- 
ward propagation of the system wave function, and the 
backwards propagation of the new Lagrange multiplier 
wave function. At this point it is worth noting that the 
presence of forward and backwards time progations is a 
general feature of the quantum kinetic theory which can 
be conveniently formulated as a propagation along the 
Keldysh-Schwinger closed-time contour [l^[l^. There- 
fore it is natural to expect that there is a connection 
between QOCT and the Keldysh contour formulation of 
the quantum dynamics. In the present work we make 
this conection explicit by re-examining the derivation of 
the expression of the gradient (or functional derivative) 
of the target functional. 

Our main simple observation is that the differentia- 
tion of a target observable with respect to a control pa- 
rameter is identical to computing a change of that ob- 
servable induced by a corresponding perturbation in the 
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Hamiltonian. Thus, the problem of calculating the gra- 
dient of the merit function reduces to a generalized form 
of linear response theory (LRT), in which the unper- 
turbed Hamiltonian is no longer static but depends on 
time. The formalism of LRT can then be directly ap- 
plied, and we straightforwardly recover the very same 
expressions that one reaches in the "traditional" way. 
However, these expressions can then be regarded as re- 
sponse functions represented by certain retarded correla- 
tion functions. We emphasize that this re-derivation is 
not a mere academic exercise, since the new interpreta- 
tion of the gradient as a response function suggests im- 
mediately the use of the known approximations to this 
object. In particular by relating the retarded response 
function to a contour-ordered correlation function we 
can apply well developed methods and approximations 
of the non-equilibrium many-body perturbation theory 
to QOCT for many-electron systems [H, [13, [HI . 

The latter is a specially important aspect, since the 
treatment of many-electron systems in notoriously diffi- 
cult; yet the direct control of electrons is an area of grow- 
ing interest, due to the advances in laser pulses of strong 
intensity and ultra-short durations, in the atto-second 
range - the scale of the electronic movements. In order to 
theoretically study a direct control of electronic motion, 
it is necessary to have a predictive (ab initio) yet com- 
putational tractable scheme, in combination with QOCT. 
Some possibilities have been recently put forward, such as 
(multi-configuration) time-dependent Hartree Fock 
and time-dependent density- functional theory Here 
we propose a new possibility, based on non-equilibrium 
many-body Green's functions theory. 

The structure of the paper is the following. In Sec. HIl 
we derive the gradient QOCT equations in the formalism 
of LRT. To make this paper self-contained, the slightly 
generalized basic LRT results needed for this purpose 
are presented in Appendix [A] Sec. IIIII elaborates on 
the equations derived in Sec. [H] by proposing a QOCT 
scheme for many-body systems, based on the Keldysh 
contour formalism and on standard approximations in 
non-equilibrium many-body Green's functions theory. 



II. THE BASIC QOCT EQUATIONS IN THE 
LINEAR RESPONSE THEORY LANGUAGE 

Let us consider a quantum system described by its den- 
sity matrix p{t) and governed, in the time interval [to, tj], 
by a von Neumann equation in the form: 

= [H[u]{t),p{t)] , 
/5(to) = Po , 
where the Hamiltonian is given by [23 |: 

H[u]{t) = n + e[u]{t)V . 



(1) 
(2) 

(3) 



a set of parameters that we will denote, collectively, u. 
The operator V represents the coupling of the system 
with an external field, e.g. if we think of an atom or 
molecule irradiated by a laser pulse, the dipole operator. 
Evidently, a particular choice of the control u leads to a 
system evolution, u p[u\(t). 

We wish to find the values of u that maximize the value 
of the expectation value of some observable A at the end 
of the propagation. In other words, we want to find the 
maximum of the function; 



G[u]^Tr{p[u]{tf)A}^ 



(4) 



In order to find the maximum, the best way is to be able 
to compute the gradient of G. The problem that we face, 
therefore, is that of finding a suitable expression for this 
gradient. 

Assuming that there is only one parameter u (the gen- 
eralization to more than one is trivial): 

BC 

— \u]= lim Au-\G[u + Au] - G[u]) . (5) 

Note that p[u] corresponds to the propagation of the 
system with the Hamiltonian given in Eq. ([3]), whereas 
p[u + Au] corresponds to the propagation of the system 
with the Hamiltonian 

de 

H[u + Au]{t)^H[u]{t) + Au—[u]V, (6) 

ou 

to first order in Au. Now we can use directly the LRT 
result introduced in appendix by making the identifi- 
cations: 

Ho{t)^HM{t), f{t)^Au^[u]{t). (7) 

ou 



Therefore, we just need to apply Eqs. (jA12p and (|A13p 
to arrive at: 



dG 



f°° de 



(8) 



where 



XA,y(^/'^) = -^e{tf^T)Tr{p{to) [i^ (ty), (t)J } (9) 

is the response function for the {A, V) operators. Inside 
the commutator, these operators appear in the Heisen- 



berg representation, defined by: 



dH{t) = U\tM)OiJ{tM) 



(10) 



The Hamiltonian piece % is static, and e[''i](t) is a time- 
dependent function whose precise form is determined by 



for any observable O, and where U{t,to) is the propaga- 
tor corresponding to the H[u]{t) Hamiltonian. Eq. ^ 
clearly manifests how the gradient is nothing else than 
a response function - albeit a generalized one. It cor- 
responds to the response of a system driven by a time- 
dependent Hamiltonian, to a modification of this Hamil- 
tonian. It remains now to see how this result is equivalent 
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to the expressions obtained in a different manner with the 
usual QOCT technique. For that purpose, we define an 
operator: 

A[u]iT)=UiT,tf)AU{tf,T), (11) 

A[u]{tf)^A, (12) 

which can also be written as the solution to the differen- 
tial equation: 



|iM(i) = -z H[u]it),A[u]{t) 
A[u]{t^^)^A. 



(13) 
(14) 



Using this new auxiliary object, and after a little manip- 
ulation of Eq. ([5]) one arrives at: 



dG 



-z I'dT ^[u]iT)Tr{p[u]iT) A[u]{t),V } 



to 



du 



(15) 

Equations ((T^ . ((Tl)) and (ITSI) . together with the orig- 
inal propagation equation for p[u]{t), are the "QOCT 
equations" , usually derived in a different way (through 
the definition of a Lagrangian function) . Algorithmically, 
the computation of the gradient is performed with two 
consecutive propagations, one forwards for the original 
system equations, and one backwards in order to obtain 
A[u](i). These propagations provide the necessary in- 
gredients to compute Eq. ([T5|) . In the next section we 
will make a link of these forward and backwards propa- 
gations to the formulation of the quantum dynamics via 
the Keldysh contour formalism. 

It is also easy to see that all variations and generaliza- 
tions of the QOCT equations naturally follow from our 
linear response approach. 



where |x[^^](0) is defined by the expression 
\x[u]{t))^U{t,tj)A\Mu]itf)). 



(19) 



Alternatively this function can be viewed as a solution 
to the following backwards propagation problem 

d 



dt 



\x{t)) ^ -^H[u]{t)\x{t)) . 



|x(i/))-i|*M(t/)) 



(20) 
(21) 



which coincides with the standard QOCT equations 
for pure states. Within the usual formalism the state 
x[i*](i)) appears as a "Lagrange multiplier" wave func- 
tion. 



B. Continuous parameters 

The case in which the control function e(i) is not pa- 
rameterized, but one does the search in the whole space 
of continuous functions, can also be treated essentially in 
the same manner. In this case, instead of a gradient we 
will obtain a functional derivative; in fact, this derivative 
is nothing else than the response function, i.e. Eq. ([5]) is 
simply: 



5G 

W) 



This can be rewritten, for the pure state case, as: 

^^=2im{xm)\vmm): 

where x[e](t) is the solution to Eqs. and (PT|) . 



(22) 



(23) 



Pure states 



C. General target functionals 



For the case of a pure state dynamics the density ma- 
trix takes the form p[u]{t) = |5'[u](t))(^[w](i)|, where the 
wave function |5'[u](t)} evolves from a given initial state 
I^M(^o)) = 1 4*0) according to the Schrodinger equation 



|l*N(0) 



-iH[u]{tmu]it)). 



(16) 



The gradient of the merit function is given by the general 
Eq. ([S]). The only difference is that now the initial density 
matrix entering the response function describes a pure 
state: po = |^'o)(5'o|- Hence Eq. ^ reduces to the form 

r) = -ie[t - t)(*o| [Anit), Vh{t)\ 1*0). (17) 

Inserting this equation into Eq. (|5]), writing the commu- 
tator explicitly, and inspecting the terms we find that the 
gradient can be written as follows 



dG_ 



2Im 



*dr|^M(T)(xM(r)|y|*M(r)). (18) 



In some cases, the function to optimize is not a simple 
expectation value of an operator A, but perhaps a more 
general expression in the form: 



G[u]^F[p[u],u]. 



(24) 



where F is a functional of the evolution of the system 
(and also perhaps explicitly of the control parameters, 
hence the second argument). Normally, this is split as: 



(25) 



i.e. the first term is the real objective, depending on 
the evolution of the system, whereas the second term 
is added in order to penalize undesired features of the 
control function, such as for example too high frequencies 
or intensities. In any case, any physically meaningful 
definition for Ji will be that in which it is a function 
of expectation values of observables. In this case the 
derivation outlined here is directly applicable, by a simple 
use of the chain rule. 
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to 



Let us reconsider the key equation for the gradient of 
the merit function, Eq. and write it exphcitly as 
follows 



FIG. 1. Keldysh contour. 



D. Time-dependent targets 

A more interesting generalization is that in which the 
function to optimize depends on the expectation value of 
the operator at all times during the propagation, and not 
only at the final time t f : Once again, this case can also be 
put in response- function language in a rather straightfor- 
ward manner. Let us consider for example the pure-state 
case: 

G[u] = fdt gitmu]{t)\Amu]it)) , (26) 

Jto 

where g{t) is some weight function. The application of 
the LRT equations leads now to: 

— M = dt drg{t) — [u]{r)xA,v{t.r). (27) 

Here the response function vi^i^) '^^ given by Eq. (1171) . 
Following the same route as in derivation of Eq. (1181) in 
Sec. niA we rewrite Eq. (^7)) as: 

— M = 2Im dr —[u]{T){x[u]{r)\V\Mu]{T)} , (28) 
where x[u](t) is defined by the following integral 

lxN(^)) = jy g{t)U{t, r)i|*M(t)) , (29) 
which can be put in the equivalent differential form: 
l-lxM(r)) = -zH[u]{rMu]ir)) - giT)Amu]iT)m 

OT 



\x[u]{tf))=0. 



(31) 



These are once again the backwards QOCT equations, in 
the case of "time-dependent targets" . 



dG 



- I pdr ^[u]ir)Tr{p{to)VHiT)AHitf)}m 



The two integrals in this equation can be composed into 
a single integral over the Keldysh contour C depicted in 
Fig. ([T|). This contour starts at to, goes froward in time to 
t f , and then comes back to the origin. Therefore by using 
the standard definition of a contour-ordered correlation 
function 



x'Xy{r,r') = -tTr{p{to)Tc AH{r)VH{T') } (33) 



where Tc is the chronological ordering operator on the 
contour C, we can cast Eq. ([32]) into the following com- 
pact form 



M = j^dr-[u]{r)x%y{tf.r), (34) 



du 



The main advantage of the representation p4p is 
that for interacting many-body systems the contour- 
ordered correlation functions can be calculated using the 
standard diagrammatic technique for non-equilibrium 
Keldysh Green's functions (see, e. g„ Refs.[ll,[23,[2ll,[2l- 
[28l ). In other words, by employing the well developed ma- 
chinery/approximations of the non-equilibrium Green's 
functions theory (NEGFT) we can express the gradients 
of the merit function as a functional of the contour or- 
dered one-particle Green's functions. 

To illustrate above statements we consider the simplest 
situation when both the control field V and the observ- 
able of interest A are represented by one-particle oper- 
ators. In this case the correlation function x'^yiifi''') 
entering Eq. (p4l) is given by: 



III. QOCT IN TERMS OF THE KELDYSH 
CONTOUR FORMALISM 

The new point of view on QOCT proposed in the previ- 
ous section naturally suggests new approximation strate- 
gies for control problems in interacting many-electron 
systems. As we will now show the QOCT equations can 
be expressed in terms of correlations functions defined 
on a Keldysh [l^l closed time contour. This allows for 
an immediate application of the powerful machinery of 
non-equilibrium Green's functions theory to the coher- 
ent control problem. 



(35) 



where K is the exact two-particle Green's function. Now 
we can take our favorite many-body approximation, such 
as Hartree-Fock, second-Born, T-matrix, random phase 
approximation (RPA), etc., to get an explicit and practi- 
cally feasible expression. For example, at the RPA/GW 
level the correlation function reduces to the two following 
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terms: 




W 




V 



T2 



(36) 



Analytically, this diagram translates to: 
xly{tf,r)=tT{AG{tf,r)VGiT,tf)} 

+ J dTidT2 j dridr2tT{AG{tf,T)h{rx)G{T, tj)} 
xWiri,T,;r2,T2)tT{h{r2)G{tf,T)AG{T,tf)} (37) 

where G(ri, T2) — G(ri, ti; r2, T2) is the one-particle con- 
tour Green's function, VF(ri, ri; r2, T2) is a dynamically 
screened Coulomb interaction, h(r) is a one-particle den- 
sity operator, and all traces are taken over a one-particle 
Hilbert space. 

Equation (1571) shows that for the practical calculation 
of the correlation function ^{tf,T), and thus the gra- 
dient of Eq. (|34| we need the countour ordered Green's 
function G and the screened interaction W. The latter 
is given by the RPA integral equation; 



(38) 

while the former is calculated by propagating the 
Kadanoff-Baym equation; [20], 



h{l) G(l, 2) = 6{l, 2) + / d3E(l, 3)G(3, 2) 



A 

(39) 

and its conjugate on the time contour. In Eq. (j39p 
/i(l) = /i(ri,Ti) is the one-particle Hamiltonian which 
also includes the Hartree potential, and the self energy is 
given by the GW diagram; 



E(1,2) = G(1,2)W^(2,1) 



(40) 



More technical details can be found, for example, in 
Ref. I27I. At this point it is worth to comment on one 
technical issue. Most currently existing implementations 
of the Kadanoff-Baym equations [26l - l28| assume that the 
dynamics starts from the thermal equilibrium state at 
some temperature T = The equilibrium initial con- 
ditions are technically convenient because they can be 
treated by a slight modification of the Keldysh contour. 



Namely, one attaches a "vertical track" going from to to 
to — if3 from the backward branch of the contour, and 
imposes antiperiodic Martin-Schwinger boundary condi- 
tions G(to — "iPjT) = —G{to,T) on the Green's function. 
If this formalism is employed then all time integrations 
in Eqs. and ([55]) are along the modified contour in- 
cluding the vertical track. However, this does not infiu- 
ence the calculation of the gradient of Eq. as it re- 
quires only the correlation function on the real time, for- 
wards and backwards branches of the contour. We would 
like to emphasize that the use of equilibrium/ground 
state initial conditions is not a fundamental restriction 
of NEGFT. It is also possible to formulate the theory for 
a general initial state [2l|, |2^ [2^ |33| although we are not 
aware of any practical implementation of this formalism. 

We conclude this section by noting the following re- 
markable fact regarding the Keldysh contour formula- 
tion of QOCT for interacting many-body systems. If the 
quantum dynamics is described within NEGFT the im- 
plementation of QOCT does not require solving any ad- 
ditional equation. All ingredients required to calculate 
the merit function gradients are already known from the 
solution of the Kadanoff-Baym equations. For example 
at the RPA/GW level of the theory one only needs to 
plug the known functions G and W into Eqs. ([55)1 and 
P4p. perform the integrations, and close the optimization 
loop. 



IV. CONCLUSIONS 

We have shown how the key equations of QOCT can 
be easily derived by employing the formalism of linear 
response theory. These equations provide the gradient 
of the target functional with respect to the external field 
which has to be optimally shaped. In the light of the lin- 
ear response interpretation, the gradient is in fact the 
response function of the driven system. First of all, 
this derivation is valuable methodologically as it explains 
the internal structure of the coherent control theory us- 
ing one of the most common techniques in theoretical 
physics, thus making QOCT more clear and accessible 
to a broad audience. In addition to that our LRT repre- 
sentation immediately suggests a reformulation of QOCT 
equations in terms of the Keldysh contour-ordered corre- 
lation functions. The theory of non-equilibrium Green's 
functions (NEGFT) may then be directly applied to de- 
rive new approximation strategies for control problem in 
interacting many-electron systems. We stress out that 
the implementation of QOCT looks especially simple, 
if the quantum dynamics is described within NEGFT, 
as it is frequently done in practice for many-body sys- 
tems. To calculate the merit function gradients there is 
no need to solve any additional equation, since all the 
required quantities are already known from the solution 
of the Kadanoff-Baym equations. Work along this line is 
in progress. 
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Appendix A: Generalized Kubo formula 

Let us consider a system governed by a total Hamilto- 
nian H{t) that is split as: 



H{t) ^ Ho{t) + f{t)V , 



(Al) 



given some real time-dependent function f{t) supported 
in the time interval [to,^/]- To formulate a generalized 
LRT we need to solve the equation of motion for the 
density matrix p{t) 

t^^m^[Ho{t)+fm,m] (A2) 

for some given initial /3(to), by considering the second 
term as a "perturbation" , while allowing the first term, 
Hq, to be time dependent. 

We search for a solution in the following form 



p{t)=Po{t)+Pi{t), 



(A3) 



where poit) solves Eq. (IA2[) with f{t) = and the initial 
condition po{to) = p{to), and pi{t) is a solution to the 
linearized equation 



^§iPi{t) = [Ho{t),Pi{t)] + [f{t)V,po{t)] (A4) 



with the initial condition pi{to) = 0. 

It is convenient to introduce a propagator U{t,t') for 
the unperturbed evolution 



(A5) 



where T is the usual time-ordering operator. Equa- 
tion (jASP is a formal solution to the equations 



i-^U{t,t') = -U{t,t')Hoit') 



with the boundary condition U {t, t) = I. 



(A6) 



Using Eqs. (|A6[) we immediately find both the unper- 
turbed density density matrix po(0 and the solution pi (t) 
of the linearized equation (jA4| : 

(A7) 



Poit) = Uit,to)p{to)U{to,t), 



Pl{t)^-t dTU{t,T)[f{T)V,po{T)]U{T,t). (A8) 



It is easy to check that pi (t) of Eq. (jA8l) is the solution 
to Eq. (jA4p . Indeed, the differentiation with respect to 
the upper limit of the r-integral in Eq. (|A8[) yields the 
second term in the right hand side in Eq. (|A4p , while the 
i-derivatives of the propagators in Eq. (jA8P produce the 
first term, [Ho{t), pi{t)]. 

Now one can calculate the change SA{t) of the expecta- 
tion value for any observable A, which is induced by the 
perturbation [the second term in the Hamiltonian ()Al|) ]: 



M(t) =Tr{pi(t)i}. (A9) 

Inserting pi (t) of Eq. (jA8P into Eq. (jA9P and rearranging 
terms we get the result 

SAit) = -^ /" dT /(T)Tr {p{to) [Anit], Vh{t)] } , 



(AlO) 

where operators 0{t) in the Heisenberg representation 
are defined as follows 

Onit) := C/(to, t)dU{t, to) = UHt, to)dU{t, to). (All) 



Equation (jA10[) suggests the definition of the (A, V) re- 
sponse function as: 

XAyit^t') = - ^')Tr {p{to) [^^(t), f//(t')] } 

(A12) 

so that: 



SA{t)^ / dTfiT)xA^{t,T). 
J to 



(A13) 



The response function of Eq. (jA13|) has the standard 
form of Kubo's formula 31]. The only minor difference 
is that for a time-dependent unperturbed Hamiltonian 
Ho{t) the Heisenberg operators, Eq. (|A11|) . are defined 
via the time-ordered exponential of Eq. (|A5I) . 

Finally, we note that the QOCT equations can also be 
derived in yet another different but equivalent manner 
by making use of the following identity for the quantum 
mechanical propagator associated to a Hamiltonian that 
depends on a parameter A: 



d 



UxitfM) 



~i I dtul{t,tf)^{t)Uxit,to) 



(A14) 

With this identity, it is straightforward to compute the 
derivative of: 

G[u] = {nu]{tf)\Amu]{tf)) 

= {'fo\Uuito,tf)\A\U^{tf,to)\^o) , (A15) 

where J7,j(i,to) is the propagator determined by the 
Hamiltonian H[u]{t). 
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